function [w,wUse] = OMP(y,X,N,err)
% OMP is the Orthogonal Matching Pursuits algorithm for regression
%
% Written by: John Paisley, Duke University, jwp4@duke.edu

magy = y'*y;
w = zeros(size(X,2),1);
wUse = zeros(size(X,2),1);
Xdot = X./repmat(sqrt(sum(X.^2,1)+eps),length(y),1);
r = y;
bool = 1;
maxIters = 50;
iter = 1;
while bool
    iter = iter+1;
    dotBas = Xdot'*r;
    [tmp,idxAdd] = max(abs(dotBas)); 
    wUse(idxAdd) = 1;
    idxLib = find(wUse == 1);
    wOpt = (X(:,idxLib)'*X(:,idxLib))\(X(:,idxLib)'*y);    
    r = y - X(:,idxLib)*wOpt;
    if iter>maxIters || (r'*r)/magy < err || sum(wOpt~=0) == N 
        bool = 0;
    end
end
w(idxLib) = wOpt;